#ifndef COMB_SNP_CAND_HPP_
#define COMB_SNP_CAND_HPP_

#include "io.hpp"
#include "util.hpp"
#include "perm_parser.hpp"
#include "sam_parser.hpp"

#include <iostream>

#include "comb_types.hpp"

/** Because of the structure of color space reads.  Not all color sequences are valid SNPs.  Some color sequences are more like errors then SNPs.  This will return true if the two color sequences are valid SNPs.
 * \param read the read sequence
 * \param ref  the reference sequence.
 */
bool two_colors(const char *read, const char *ref);

/** Decide if the read contains a valid snp block of length <= len.
 * \param read the read in color space.
 * \param ref the reference in color space.
 * \param len the length to be tested.
 */
bool legal(const char *read, const char *ref, unsigned int len);

/** Use the genome to add sequence information to snp blocks. 
 * \param seqs the genome
 * \param snp_cands the SNP candidates.
 */
void attach_seqs_to_snp_cands(std::map<std::string, std::string> &seqs, cand_vector_t &snp_cands);

/** Calculate the number of reads that map to each location, and mismatches at each location.
 * \param options Options passed to comb
 * \param subs returned array of substitution counts for each position.
 * \param totl returned array of total mapping to each position
 * \param ref_seq_len genome sizes.
 */
void find_subs_totl(prog_options &options, map_str_vafl_t &subs, map_str_vafl_t &totl, map_str_ui_t &ref_seq_len);

/** Using the subs/totl arrays calculate SNP candidate sites.
 * \param options options given to comb
 * \param cand_blocks SNP candidates
 * \param snps SNP candidates
 * \param subs the substitution counts per base
 * \param totl the total mapping count per base.
 */
void find_snp_cands(prog_options &options, cand_blocks_t &cand_blocks, cand_vector_t &snps, map_str_vafl_t &subs, map_str_vafl_t &totl);

/** Calculate the number of reads that map to each location, and mismatches at each location. In sequence space.
 * \param options Options passed to comb
 * \param subs returned array of substitution counts for each position.
 * \param totl returned array of total mapping to each position
 * \param ref_seq_len genome sizes.
 */
void find_seq_subs_totl(prog_options &options, map_str_vafl_t &subs, map_str_vafl_t &totl, map_str_ui_t &ref_seq_len);

/** Using the subs/totl arrays calculate SNP candidate sites.
 * \param options options given to comb
 * \param cand_blocks SNP candidates
 * \param snps SNP candidates
 * \param subs the substitution counts per base
 * \param totl the total mapping count per base.
 */
void find_seq_snp_cands(prog_options &options, cand_blocks_t &cand_blocks, cand_vector_t &snps, map_str_vafl_t &subs, map_str_vafl_t &totl);

/** Driver for "comb mapcand"
 */
void do_mapcand(int argc, char **argv);

/** Driver for "comb mapcount"
 */
void do_mapcount(int argc, char **argv);


/** Driver for "comb snpcand"
 */
void do_snpcand(int argc, char **argv);
#endif
